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Abstract. The first, self-consistent calculations are presented of the cosmological. 
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fvq i H2-dissociating UV background produced during the epoch of reionization by the 

I^^V I sources of reionization. Large-scale radiative transfer simulations of reionization trace 

f^*) ' the impact of all the ionizing starlight on the IGM from all the sources in our simulation 

QO , volume down to dwarf galaxies of mass ^ 10® Mq, identified by very high- resolution 

N-body simulations, including the self-regulating effect of IGM photoheating on dwarf 
galaxy formation. The UV continuum emitted below 13.6 eV by each source is then 
transferred through the same IGM, attenuated by atomic H Lyman series resonance 



^ ' lines, to predict the evolution of the inhomogeneous radiation background in the 



Lyman- Werner bands of II2 between 11 and 13.6 eV. On average, the intensity of this 
Lyman- Werner background is found to rise to the threshold level at which dissociation 
suppresses H2 cooling and star formation inside minihalos, long before reionization is 
complete. Spatial variations in the Lyman- Werner background are found which result 
from the clustering of sources associated with large-scale structure formation, such 
that intensity fluctuations correlate with matter density fluctuations. As a result, the 
Lyman- Werner background rises to the threshold level for II2 suppression earlier in the 
vicinity of the reionization sources and their H II regions. 



1. Introduction 

Simulations suggest that the first stars in the Cold Dark Matter ("CDM") universe 
formed inside minihalos of mass M ~ 10^~^Mq at 2 > 20, when H2 molecules 
cooled the primordial, metal-free halo gas and gravitational collapse ensued (e.g. 
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Bromm and Larsonll2004j and references therein) . This critical role of H2 molecules as the 
primary coolants responsible for triggering the gravitational collapse that caused stars 
to form inside minihalos was limited, however, by the fact that H2 can be dissociated 
by absorbing UV radiation in the H2 Lyman- Werner ( "LW" ) bands in the energy range 
11 — 13.6 eV. In the presence of a high enough LW-band radiation intensity, Jlw, the 
H2 abundance would have been too low to cool the gas sufficiently to form stars. The 
threshold level of the intensity, (Jlw) threshold, above which rainiha lo star formation was 



suppressed is still uncertain. Early estimates ( iHaiman et al.ll2000l : henceforth, "HAR") 
found that (Jlw) threshold depended upon minihalo mass and redshift. In terms of the 
dimensionless quantity, Jlw, 21 = Jlw/(10~^^ ergs~^ cm~^Hz"^ sr~^), they found that, 
to suppress H2 for all minihalo masses, (Jlw, 21) threshold ~ ^^~^ ^'^ ^ ^^^ required, as 
redshift varied from 2; ~ 10 to 50, respectively. Later estimates, including those based 
upon 3D, numerical gas dynamical simulation of minihalos evolving in the presence of 
a LW background from more r ealistic cosmologica l initial conditions, found a similar 



range of threshold values (e.g. iRicotti et al.ll2002l ; lYoshida et al.l l2003l ; lYoshida et al 



2003). 

Once stars began to form inside minihalos and also inside the rarer, more massive 
halos with M > 10^ Mq virial temperatures above lO^K, in which radiative cooling 
by atomic hydrogen was possible even without H2 (to cool the gas down to lO'^K, at 
least, and initiate gravitational collapse), a rising diffuse LW background would have 
been inevitable, however. Starlight at energies below 13.6 eV, the ionization potential of 
atomic hydrogen, would have been largely free to escape from these star forming halos 
into the intergalactic medium ("IGM"), while the ionizing radiation above 13.6 eV was 
partially absorbed by the neutral hydrogen within the halos and thereby reduced by 
an escape fraction, /esc- In that case, given the value of this /esc and the ratio of the 
number of ionizing photons to the number of H2-dissociating photons released by the 
stars, A^i/A^LW, which depends on the mass function and spectra of the stars, the rise of 
the cosmic LW background can, in principle, be related to the rise of the diffuse ionizing 
background. The latter is believed to have been responsible for the reionization of the 
intergalactic medium completed by redshift z >6. Estimates by HAR showed that the 
sources of this reionization would have caused the mean LW intensity in the IGM to 
exceed (Jlw, 21) threshold ^°^S before reionization was complete. 

This outcome is expected on quite general grounds, in fact. The mean number of 
LW photons in the background per H atom is related to the LW intensity according to: 
/nLw\ /47r /•i=^-^<=^ J, 



V nil 



( — / -r^du ) /nn, 

V c Jn.5ev hu J 

^12:6^^'^^^/^«' 
1.05 X 10-^ Jlw, 21 /nn, 

-3. /1 + ^ 



~6.3xlO-VLw,2i(^pl . (1) 

Hence, if (Jlw, 21) threshold < 1' this implies that (wlw/'^h) threshold < 6x 10"^(l+2;)2i^ < 1. 
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The mean number of photons in the LW background per H atom at a given epoch is 
determined by their rate of emission per atom, integrated over time, reduced by a factor 
which accounts for attenuation and redshifting after emission. As discussed in § El 
photons emitted in the LW bands below 13.6 eV are removed from the background 
when they redshift to the frequency of the nearest H atom Lyman series transition from 
level ra > 3 to n = 1. Only those sources within the distance from which photons 
emitted at the Ly7 frequency would be received at the Ly/3 frequency can contribute 
at all to the LW background at a given point. On average, roughly a third of the LW 
photons emitted within this horizon will survive their tripj. The same sources also emit 
UV photons with energies above 13.6 eV which are destroyed en-route by photoionizing 
H atoms. Reionization required that at least one ionizing photon was released into the 
IGM per H atom by the end of reionization. Let ^lw and ^i be the total number of LW 
and ionizing photons, respectively, released per H atom into the IGM up to some time. 
The mean number of photons per H atom in the LW background at that time is then 
roughly given by 

'^H 3 V 6 / 

For the stars believed to be responsible for reionization, the intrinsic ratio of ionizing 
UV photons to LW photons released, A^j/A^lw; ranged from Aj/Alw ~ 15 (Pop III stars, 
high mass) to Ai/Alw ~ 1 (Pop II stars, Salpeter IMF), while the escape fraction of 
ionizing photons was some /esc <^ 1. In that case, we can write 

and, therefore, if (Jlw, 21) threshold < 1, equations ([H) and ([3]) imply that (wlw/'^h) reaches 
the threshold level when ,^i is only 

I^OlW, threshold — ^ I — I I ITT I /esc ^ 1- (4j 

V ^H /threshold V^^LW/ 

Accordingly, since reionization was not complete until the condition ^i > 1 was reached, 
the LW threshold for suppressing H2 in the minihalos must have been reached long before 
the end of reionization. As a result, as HAR suggested, the minihalos were generally 
"sterilized" before they could contribute significantly to reionization. 

What level of LW background is required to suppress minihalo star formation 
when the minihalos are also directly exposed to ionizing radiation, as well, is a 
more complicated question to answer. If a minihalo forms in a region of the IGM 
which is already photoionized, the pressure of the photoheated IGM there prevents it 
from collapsing gravitationally into the dark-matter dominated halo. Those minihalos 
are missing their baryonic componen t, therefore. This phenomenon, sometime s 



referred to as "Jeans-mass filtering" (jShapiro et al.l . ll994J : iGnedin and Hull . Il998l ). 



J This factor 1/3 comes from the survival rate from the attenuation by Lyman series resonace hues 
given approximately by Jq^^ drosfmodifos)/ Jq^^ dros, where the terms used are described in § 12.2 
Compare this expression to equation ([29|) and also see Figure [H 
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ensures that H2 cooling and star formation do not occur inside niinihalos in the 
ionized regions of the IGM unless those minihalos had formed there prior to the 
arrival of the ionizing radiation. The impact of both LW and ionizing radiation 
(including X-rays) on pre-existing minihalos insi de H II regions i s a sub j ect of ongoing 
work, beyond the 
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20081 ). The presence of dissociating radiation always implies 



the potential for limiting the H2 abundance and, with it, the cooling required to make 
stars form. This is true even for the atomic-cooling halos above the minihalo mass 
range, although the level required for suppression might be higher. We shall focus here, 
however, on the rise of the LW background and its spatial variations contributed by the 
dominant sources of reionization, but leave the question of how the intensity impacts 
star formation for future studies. 

Previous estimates of the cosmic LW background were limited to the mean 
background and were based upon a homogeneous approximation. These calculations 
assumed that the sources and the IGM were uniformly distributed, with uniform 
emissivity, given either by analytical approximation (HAR) or by summing over the 
sources found in small-box simulations, too small to account for the large - scale clustering 



Ricotti et al. 



2002 



Yoshida et al. 



of so urces or to follow global reionization (e.g 
20031). Background of L ya p umping radiation was co n sidere d semi-analytically by 
Barkana and Loebl (120051 ) and iPritchard and Furlanettd (120061 ) . including the effects 
of fluctuations in the background in a linear approximation. The focus of these papers 
is on the Lya intensity originating from photons that are absorbed by hydrogen Lyman 
resonance lines and converted into Lya photons, which then radiatively mixes the 21 
cm levels to drive the spin temperature to the gas kinetic temperature (Wouthuysen- 
Field effect), while we are interested in photons that remain unattenuated by Lyman 
resonance lines and affect H2 abundance through photo-dissociation. Note that the 
horizon for sources responsible for the fluctuating Lya background is considerably larger 
than that responsible for the LW background. 

Here we present the first self-consistent radiative transfer calculations of the 
inhomogeneous LW background produced by the same sources which reionized the 
universe in a large-scale radiative transfer simulation of reionization. This problem 
presents a formidable computational challenge. The horizon for seeing LW photons is 
~ 100 comoving Mpc, much larger than the size of typical H II regions (~ 10 Mpc). 
Since the mean free path for LW photons (~ 100 Mpc) is much larger than that for H- 
ionizing photons we must account for sources distributed over large volume and lookback 
time. Finally, the LW band photons are attenuated as they redshift into H-atom Lyman 
series resonance lines as they travel across the IGM. Hence, it is necessary to perform 
a multi-frequency radiative transfer calculation from each of the millions of sources 
in a cosmological volume larger than ~ (100 Mpc)^, integrated along the light cones 
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from each source to every observation point they intersect, which is computationally 
prohibitive. As we shall show, a good approximation is possible which reduces the multi- 
frequency calculation to an equivalent gray opacity calculation. As a result, we are not 
only able to derive the evolution of the rising globally-averaged mean LW background 
during the epoch of reionization (EOR), but also to map its pattern of spatial variations 
over time. In § [2l we describe how continuum radiation emitted in the LW range below 
13.6 eV is transferred through the IGM along the light cones from sources to observers, 
and how we solve this problem numerically. In § [3], we apply this method to one of our 
rec ent large-scale radiative transfer s imula tions of self-regulated reionization, described 



m 



Iliev. Mellema. Shapiro and Pen] (120071 ). We compare the mean LW background 
evolution thus derived numerically with the homogeneous universe approximation and 
describe the spatial fluctuations of the LW background in some detail. Our conclusions 
are summarized in § HI Some of our results on the in homogeneous LW background 



during the EOR described here were first summarized in lAhn et al.l (120081 ). 

2. Radiative Transfer of the LW Background 

2.1. Basic Equations 

Let us first briefly describe how the inhomogeneous LW background can be 
calculated. Consider radiation sources distributed inhomogeneously. The mean intensity 
Jui^ohs, -^obs, t'obs) at obscrved frequency i^obs at some comoving position Xobs at redshift 
2;obs is given by 

'-'i/l.^obs; 2^obs; '^obsj "j / ^ -f u,s (.Xpbs , ^obs ; ^obs ) ; V<3J 

s 

where F^^^ is the flux received at (xobs, ^obs, t'obs) that was emitted at (xg, Zg, ^^s) by a 
source (denoted by subscript s), where 

^s ^ 1 + ^s ,„v 

^obs -L ~r •Z^obs 

The position and redshift, (xs,-2;s), of a source are related to those of the observer at 
(xobs, ^^obs) by the fact that the signal emitted at the epoch Zs must reach the position 
Xobs at the epoch Zohs, travelling at the speed of light while the universe expands. We 
express this implicitly by writing the comoving separation, Tos, of the source and observer 
as follows: 

/"*^'°''^ cdt _ r= dz 

Jt{z.) ait) 4^^ Hiz) 

If we specialize to the case of interest here [z > 6) in which the Hubble parameter 
H{z) is given by the high-redshift limit for a flat universe with cosmological constant, 
equation l^^ can be integrated to yield 

ros = 2cH,'Q;^'/'[{1 + zobs)-'/' - (1 + z,)-'/% 

-1/2 



'^obs 



jy^ 



(8) 
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using equation 1^. 

The differential flux, F^,,,, received at (xobs, ^^obs, ^'obs) from a source of differential 
luminosity L^, emitted at (xg, Zg, z/g) is given by 

IP / N L^{u = u,) f 1+Zg \ 

47riJ£(2:obs, Zs) \l + Zohs J 
where the factor {1 + Zg) / {1 + Zo^^g) reflects the fact that the observer sees the differential 
frequency interval reduced by redshift relative to the emitted interval. Here D^^Zohs, Zg) 
is the luminosity distance given by 

^'(--' - fe) (iT^) • (^°' 

where the factor (1 + Zg)/{1 + Zohs) takes account of the redshifting of photon energies 
and arrival rates, which reduce the observed flux by [(1 + 2;obs)/(l + -^^s)]^- The optical 
depth Tj^^j^^ depends on the observed frequency z/obs according to 

Tu,i,,= pb{^,z)K^,{yi,z)cdt, (11) 

where v' = i ^J" j Vohsi pb{^i z) is the baryon density at (x, 2;), K^i{yi, z) is the opacity 
at (x, z) to photons of frequency v' , and where x and z are the position and redshift 
of photons travelling along the line of sight which were emitted at (xg, Zg) and will be 
received at (xobs, ^^obs)- 

The optical depth of the IGM to continuum UV photons in the LW range between 
11.2 eV and 13.6 eV is predominantly due to resonant absorption by neutral H atoms in 
the Lyman series lines with upper states n = z, i > 3 (HAR). The mean optical depths 
in these lines can be written as follows: 

osc,i a^^p^^ ^ / .^^ I j —^ ] TGP.a -^ 1.8/osc,i^GP,a for 1 > 1, (12) 



where /osc,i and /osc,a are the oscillator strengths and Ui and Va are the frequencies, 
for lines of upper states n = i and 2, respectively, and where rGp,a is the familiar 
Gunn-Peterson optical depth in the Lyman-a transition: 



V0.044y Vo.ry vo-27y V 21 y ^ ' 

where tt-hi is the mean number density of neutral hydrogen and x-m is the mean neutral 
fraction of the IGM. Since we are most interested here in the early phases of the 
EOR, xhi = 1 — Xh tt ~ 1. As such, rop_n. ^ 10^ ^ 1. Even inside H II regions. 



however, xhi ^ 10 (llliev. Shapiro. McDonald. Mellema and Pen! . 120071 ). so rGp,Q ^ 1 



in general. As i increases, fnsr.i decreases (e.g. fnsr.i = 0.079 , 0.029, and 0.014 for z = 3, 



4, and 5, respectively (e.g. iHohne and Zimmermannl . Il982l ). so there is some Vax such 



that Tj < 1 for i > Vax- For example, if mean IGM density is assumed in the ionized 
region with xri ^ 10-^ Vax = 8, 7, and 6 with /osc = 0.0032, 0.0048, and 0.0078 at 
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z = 20, 15, and 10, respectively. Note, however, that cosmic reionization occurs in 
an inside-out fashion, suc h that overdense regions are ionized ea rher than the mean or 



underdense regions (e.g. Ilhev. Mellema. Shapiro and Peru 120071 ). The effective imax in 



ionized regions, therefore, can be much larger than the estimates above. According to 
HAR, Zmax ^ 150 in the neutral IGM. It is a good aproximation, therefore, to assume 
that the Lyman lines are optically thick for all i, since the frequency range over which 
the lines are not optically thick is a negligible fraction of the LW range below 13.6 eV. 
The opacity of the IGM due to LW band absorption by II2 is relatively unimportant 
by comparison, since the H2 concentration in the IGM is small even before the LW 
background rises to suppress it. The pre-reionization H2 concentration is ~ 10~^ 



( IShapiro et al.l . Il994j ). and so tlw < 1- We shall, henceforth, neglect this source of 
opacity. We refer the reader to ^for more rigorous justification. 

The expected number of computational operations required to evaluate 
equations ([5]) and ([9]) is A^s ■ A^g ■ iVf, where As is the number of sources, Ag is the 
number of grid cells on which Jj, is calculated, and A^f is the number of frequency-bins 
which are required to resolve the frequency dependence of r,^^^^ adequately. While these 
equations are straightforward, evaluating them numerically in a brute-force way can be 
prohibitively expensive in computational resources. For example, the H Lyman series 
resonance lines are optically thick to photons in the LW bands frequency range, and one 
should consider A^f ^ 1 to account for this effect properly. Currently, a full 3D, multi- 
frequency radiative transfer calculation is not feasible for the problem of interest. The 
effective number of sources Ag in our problem can be as large as 10^ due to the large size 
of the LW horizon, and we need about A'g > 10^ grid-points to produce a statistically 
significant result. Cosmic reionization simulations, by contrast, do not require a multi- 
frequency operation, and yet, these simulations have only just become feasible recently 
with the help of massively-parallel computers. 

The following sections describe how we overcome this technical difficulty in 
calculating the LW background by reducing A^f to 1, even though the net result becomes 
equivalent to a full multi-frequency radiative transfer calculation. We further describe 
in detail how we sum individual -F,^,s's, taking full account of the effect of redshifting 
and of the finite light-crossing time between sources and observers. 

2.2. Attenuation of B.2 Dissociating Photons from a Single Source: the "Picket-Fence" 
Modulation Factor 

For an inhomogeneous distribution of sources, we need to calculate the attenuation 
of continuum photons emitted in the LW energy range 11.2 - 13.6 eV separately for 
each individual source, by hydrogen Lyman line resonance absorption and subsequent 
cascades along the light cones from the source to every observer. Consider a source 
emitting continuum radiation at frequency Ug at redshift z = Zs- As the photon travels 
toward the observer, it is absorbed when its frequency redshifts into an H Lyman 
series resonance line and, some of the time, the decay of the excited state replaces 
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the original photon with photons at frequencies below the range of the LW bands. If 
the original photon is resonantly scattered, it is quickly reabsorbed, until all resonant 
photon s eventually tu r n int o low-frequency photons below the LW bands. For this 



reason, 



Haiman et al.l (119971 . HRL, hereafter) and HAR assumed that, whenever the 
photons emitted in this range redshifted into one of the H Lyman resonance lines, 
they were completely attenuated and turned into low-frequency photons — mostly Lya 
photons — out of the LW range. From the observer's viewpoint, this leads to a series of 
"dark screens" , defined as sharp boundaries beyond which the observer cannot see any 
sources contributing LW intensity. These are marked by the maximum redshifts z^^^^ 
defined by 

J- + ^max ^j / 1 /I ^ 

J- ~r ■Z^obs ^obs 

where i>i is the frequency of the Lyman line closest to the observed frequency t'obs from 
above. For a homogeneous universe with spatially-uniform emissivity, this results in the 
well-known "sawtooth" modulation of a uniform, isotropic LW background spectrum 
(HRL; HAR), an example of which we have plotted in Figure [T] for a homogeneous 
ACDM universe with flat-spectrum sources. 

We shall make the same assumption here, that all LW photons are completely 
attenuated once they redshift into an H Lyman resonance line with upper state n > 3. 
However, we cannot limit ourselves to the homogeneous universe approximation. Since 
our objective is to consider contributions from individual sources that are distributed 
inhomogeneously, we must, instead, calculate how continuum photons emitted by each 
source will be attenuated by hydrogen atoms which they encounter along the particular 
line of sight that connects them with a given point of observation. We shall describe the 
attenuation of an individual source here in what follows. We shall then describe in § 12.31 
how we sum over these individual source contributions to obtain the spatially-varying 
LW background intensity. 

For the homogeneous universe in which the observed spectrum is transformed by the 
sawtooth modulation shown in Figure [1], this spectrum is the result of superposing the 
spectra of sources distributed continuously in lookback time along the line of sight. In 
that case, a given observed frequency combines the effect of photons emitted at different 
lookback times which experience different amounts of redshifting before reaching the 
observer. It is natural, then, to describe the modulation from the viewpoint of the 
observer, with the edge of each saw-tooth corresponding to the "dark screens" at the 
frequencies of the H Lyman lines, and the spectrum to the red side of a given line 
originating in the past from sources nearer than the distance to the corresponding 
screen. When we consider, instead, the spectrum of a single source in an inhomogeneous 
universe, there are also "dark screens" beyond which LW photons cannot pass, but these 
screens are best described from the point of view of the source. In that case, as a photon 
travels from the source, it will survive only until it encounters a "dark screen" in its 
future, as it redshifts into the nearest Lyman line. From the viewpoint of a source, 
the dark screens are located at the maximum radii that photons emitted at different 
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Figure 1. The observed "sawtooth" modulation of the uniform, isotropic radiation 
background observed in the UV range of the LW band range for a homogeneous ACDM 
universe with flat-spectrum sources, caused by H Lyman hue opacity of the IGM (see 
text), for redshifts z = 19.2 (bottom hues, black), z = 15.7 (middle lines, blue) and 
z = 9.9 (top lines, red). The horizontal (dashed, corresponding colors) lines show the 
unattenuated mean intensity levels at these redshifts. A spatially uniform emissivity 
is assumed which evolves in time in proportion to the collapsed fraction of the matter 
density in halos massive enough to be sources of reionization, as described in § |3l Also 
plotted in vertical lines are the locations of relevant LW bands, when H2 are assumed 
to be in the ground electronic state A"^E+ with v" = and J" = 0, 1. The height 
of these lines corresponds to log(0.01 x fosc), where f nsr. is the oscillator s treng th of 
Lyman (black) and Werner (green) bands compiled bv lAbgrall and Rouefj ( 1989 ). 



frequencies can travel from a source. Instead of equation (fT4l) . these radii are marked 
by minimum redshifts 2;min defined by 



(15) 



l + ^s 

where i>i is the (observed) frequency of the Lyman hne closest to the emitted frequency 
z/g from below. A screen corresponding to z/j has the radius rj(z/s), where rj(z/s) = 
ros{j^i', J^s, Zs). This introduces finite frequency gaps in the transmitted spectrum, which 
as a consequence resembles a "picket fence", with the intensity unattenuated between 
the dark gaps, as follows. 

Consider a dark screen located at comoving distance rj(z/j+i) from the source defined 
by 



ri{l^i+l) = Tosil'U l^i+l, Zs), 



(16) 



where a photon emitted at z/j+i from a source is redshifted into z/j. At this location, 
all the photons with z/g > z/j are completely attenuated in the following way. First, a 
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photon with z/g = z/j+i will be redshifted into Ui and attenuated at rj(i/j+i)l§|. Any photons 
with z/g in between Ui and z/j+i will then be redshifted into Ui and attenuated at some 
distance shorter than ri(z/j+i). Since no photon can cross a Lyman line frequency as it 
is redshifted, the observed spectrum will be completely black inside a trough between 
Ui and z/j+i. For z/j+i < z/obs < ^i+2, this will happen at rj+i(z/j+2)- Because 

ri+i(z/i+2) < ^^il^^i+i), (17) 

the observed spectrum at rj(z/j_(.i) will also have a trough from z/j+i to z/j+2- This way, 
all the photons with z/obs > '^i are completely attenuated at Tos = rj(z/j+i). On the other 
hand, because rj_i(z/j) > rj(z/j+i), a photon with z/g = z/j has not fully redshifted into 
z/j_i but only into z/obs (n(«^i+i); t^i, ^s) at ri(z/i+i), where 



'/obsV OSl ^S! -^sj — "s 



1 + 



2c 



;i + ^s) 



1/2, 



(18) 



obtained from equation ([8]). Photons with the observed frequency ranging from z/j_i 
to i^ohs{^i{^i+i)', ^i, Zs) will have reached rj(z/j+i) without attenuation. Therefore, the 
spectrum will have full transmission for z/j_i < z/obs < t'obs('"i('^i+i); '^i, 2;s), while there 
is a completely black trough for z/obs ('"j('^i+i); '^i; -^^s) < t'obs < '^i- Similarly, the next 
lower-energy interval, defined by z/j_2 < t'obs < ^«-i, will have full transmission for 



T^i-2 < z^obs < i^ohs{ri{ui+i); z/i_i, 2;s) and a trough for z/obs (n(z/i+i) 



-1) 



^s) < i^obs < 



z/j_i. As a result, the observed relative flux is affected by what we call the "picket- 
fence" modulation, as depicted in Figure El 
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Figure 2. Transmitted spectrum from a single source: picket-fence modulation, {left) 
Relative flux observed at r^ii^e): where a photon emitted at Ly5 (Lye) frequency is 
redshifted into Ly4 (LyS) frequency, (right) Relative flux observed at rio(i^ii), where 
a photon emitted at Lyll frequency is redshifted into LylO frequency. We use a 
convention that the energy of Lyz photon is 13.6 eV (l — l/i^)- The loocation of 
several Lyman resonance lines Lyi is shown in vertical lines and denoted by i in both 
panels. 



Note that, at a fixed distance, z/obs oc z/^ from equation flTSl) . A gap appearing 
between observed frequencies Uj and z/^+i at some Tos is also proportional to z/j+i. If we 

§ A more accurate description is that a photon with frequency slightly smaller than Vi+i is redshifted 
into Vi and attenuated at distance slightly shorter than ri(i'i+i), because a photon with Vs = t^i+i is 
attenuated on-site right after being emitted. However, we give the description in the text for simplicity. 
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call the size of this gap Ai/gap, j+i, we then have 

^^gap,j+l ^ ^j+l 

as long as the range of frequency, [uj, I'j+i], is not fully covered by Az/gap,j+i. 
We define the picket-fence modulation factor /mod by 



(19) 



/ll.SeV d{hUohs) 



/.od(ro. ^s) - (exp(-r.„J) = '''■''''' ^SZ, ""^ (20) 



which is a function only of the comoving distance Tos and the source redshift Zg. We 
choose 11.5 eV as the niinimum energy of interest, because the dissociation rate for the 
LW bands at hu < 11.5 eV is negligible compared to that for those at hu > 11.5 eV 
(see, e.g.. Figure 1 of HAR). When H2 molecules absorb photons in the LW bands, 
in general only about 15% of these excitations lead to dissociation of H2. When a 
single source is observed at some comoving distance Tos, some LW bands will be excited 
by fully transmitted photons, which results in dissociation about 15% of the time, 
therefore, while other bands will not, because they reside in a trough. For sources 
of interest here, we can approximate L,, by a flat spectrum whose amplitude is given 
by the frequency-averaged luminosity (L,^) = /^j^gj^y (i(/iz/s) Lj,(z/s)/(2.1eV). This is 
a fairly good approximation, unless the spectrum is unusually steep in the narrow 
energy range of interest, [11.5 — 13.6] eV. For example, both black-body spectra with 
T = [5 — 10] X 10^ K and a power-law spectrum with L^, oc z/~^ have a maximum deviation 
from (Ljy) which is smaller than 10% in this energy range. In that case, because the 
LW bands are almost uniformly distributed in frequency, the true dissociation rate will 
be almost identical to that obtained by assuming that all the H2 LW lines experience 
the frequency-averaged value of the LW intensity after the picket-fence modulation. 
Therefore, in addition to the geometrical dilution of the incident flux, the H2 dissociation 
rate will be suppressed in proportion to /mod- This /mod is just the fraction of the total 
frequency interval from 11.5 to 13.6 eV observed at Zohs from a source at z^ at comoving 
distance Tos occupied by the full transmission windows in between the dark troughs, as 
described above. Hence, 

j ^ ' 

The picket-fence modulation factor is a key ingredient in alleviating computational 
difficulties which would have arisen due to a multi-frequency calculation. We have 
calculated /mod numerically and found a simple fitting formula which fits the true values 
within a 2% error (see Figure [3]): 

_ r 1.7 exp [-(reMpc/116.29a)°-''] - 0.7 if reMpc/« < 97.39 

/mod - \ ^ -c I ^^ nr. ^^^> 

\ if TcMpc/a > 97.39 

where rcMpc is Tos in units of comoving Mpc, and a is a scaling factor given by 



•^-^TTi oi Hrl ■ P^' 
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We call ^3(1/4), which is equal to 97.39 a cMpc, the "LW horizon" tlw- This is the 
maximum comoving distance from a source that an H2 dissociating photon can reach, 
corresponding to the distance from which the redshift produces the maximum frequency 
difference possible between two adjacent lines in the Lyman series (as long as we restrict 
the observed energy range to [11.5 - 13.6] eV). Note that all the dark screen distances 
are scaled by a, which increases as Zs decreases. 




100 



Figure 3. Picket-fence modulation factor /mod as a function of comoving distance 
''cMpc in units of Mpc. True values at selected radii r.; (1/^+1) (open circle) and a fit 
(solid curve) are plotted, a, given by equation (f23|) . is a distance scaling factor which 
depends on redshift. When rcMpc = 97.39a, /mod — 0, which sets the "LW horizon" 
for H2 dissociating radiation from a source. 



2. 3. Intensity of H2 dissociating photons from multiple sources 

The path of a photon in the expanding universe follows a null geodesic. The Friedmann- 
Robertson- Walker metric for a homogeneous, isotropic universe is given by 

ds^ = dt" - a" it) [dr^ + r^d^^] = a^{t) [dr^ - {dr^ + r^d^^)] , (24) 

where we adopt natural units with c = 1. In that case, the comoving distance travelled 
by a photon since its emission is given by setting ds'^ = in equation (12^ . solving for 
dr [dVt = 0), and integrating over time to yield the conformal time r, defined by 



dr 



dt/a{t). 



(25) 



We use this fact to construct the world lines of sources and their radiation, showin 
in Figure HI If our choice of space-time coordinates is the comoving distance and the 
conformal time, then null geodesies make straight lines at a 45° angle. World lines of 
the sources, on the other hand, will be close to straight lines, parallel to the conformal 
time axis. For simplicity, we will neglect the small peculiar motion of halos. 

We must account for the finite light-crossing time for light from sources to reach 
an observer, because these are distributed over a truly cosmological volume and the 
population of sources can vary significantly over the lookback time corresponding to 
tlw ~ lOOcMpc, due to the rapid evolution of cosmological structure. The conformal 
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space-time diagram of sources mentioned above becomes a useful tool for this task. At a 
given redshift we draw past light cones from an observing point, which have a maximum 
length equal to the LW horizon length, tlw- When the world line of a source intersects 
one of these past light cones, we add its flux contribution to the mean intensity at the 
corresponding observing point (see Figure H]). The fact that At = Tos, where At is 
the conformal lookback time to a source at comoving distance r^s from the observing 
point, makes it easy to find these intersecting points as well as the source redshift. The 
conformal time interval, tlw, which corresponds to the LW horizon, tlw, determines 
how far back in look-back time a given observer cell at a given epoch 2;obs niust extend 
its past light-cone to look for contributing sources. Accordingly, this operation requires 
that the past light-cones extend back through a number, risteps? of time steps. At, equal 
to Tun/ At. 

After a contributing source is found, its frequency- averaged LW flux observed at 
the given Xobs and Zobs is evaluated using equation ([9]), replacing L^, by its average over 
the LW band frequency and replacing exp(— Tobs) by /mod using equations (l22l) and (l23l) . 
We sum fluxes from all the sources (denoted by the subscript s) observed within the LW 
horizon. 

3. The Inhomogeneous LW Background from a Simulation of Cosmic 
Reionization 

3.1. Illustrative Case: Self -Regulated Reionization 

As an example, we apply the methodology for calculating the fluctuating LW back- 
ground described in the previous sections to one of our large-scale N-body and radiative 



transf er simulations of cosmic reionization presented in llliev. Mellema. Shapiro and Pen 



(I2OO7I . henceforth "IMSP"). The cosmological structure for mation and eyoluti on is fol- 



lowed with a particle- mesh N-body code called PMFAST (JMerz et al.ll2005l ). These 
N-body results then provide the evolving density field of the IGM (coarsened to a 
lower resolution than the original particle-mesh grid in order to make the radiative 
transfer feasible) and the location and mass of all the halo sources, as input to a sep- 
arate radiative transfer simulation of inhomogeneous reionization. The latter simula- 
tion is performed by our C^— Ray (Conservative, Causal Ray- Tracing) code, a grid- 
based, ray-trac i ng, ra diative transfer and nonequilibrium chemistry code, described in 



Mellema et al.l (120061 ). The ionizing radiation is ray-traced from every source to every 
grid cell at a given timestep using a method of short characteristics. The code is ex- 
plicitly photon-conserving in both space and time, which ensures an accurate tracking 
of ionization fronts, independent of the spatial and time resolution, even for grid cells 
which are optically thick to ionizing photons and time steps long compared to the ion- 
ization time of the atoms, with correspon dingly great gains in efficiency. The code has 



been tested against analytical solutions (jMellema et al.l l2006l ) and, in direct compari- 



son with other radiative transfer methods, on a standardized set of benchmark problems 
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comoving distance (Mpc) 

Figure 4. Conformal space-time diagram of radiation sources and the past light cone 
of an observer, used to identify which of the sources in the N-body simulation volume 
in the past emitted light just now reaching the observer at a given time. Radiation 
sources are created discretely in time in the N-body simulation results — i.e. source 
catalogue is constructed at each output time (dotted horizontal lines). The location 
of each source is assumed to stay constant during each time step (shown as solid, 
dotted, and short-dashed vertical lines) as are the source luminosities. Photons follow 
null geodesies truncated at the LW horizon 7'lw, and an observer at point P will see 
sources whose world lines intersect with the past light-cone. The flux contributed by 
these sources are determined by where these intersecting points lie along the time axis. 
We use a coordinate system composed of comoving distance and conformal time, for 
computational ease. For example, the conformal lookback time Ar (Mpc) from point 
P to a source at point A, which determines the emitting redshift and the source flux, 
is easily obtained once the comoving distance r (Mpc) to the source is known, because 
Ar = r. 



(Illiev. Ciardi. Alvarez. Maselli. Ferrara. Gnedin. Mellema. Nakamoto. Norman. Razoumov. Rijkhorst. Rit 



20061). 



We simulated the ACDM universe with 1624^ dark matter particles of mass 10^ M©, 
in a comoving simulation volume of (35 h~^ Mpc)^. This allowed us to resolve (with 100 
particles or more per halo) all halos with mass of 10^ Mq or above. This is roughly the 
minimum mass of halos which can radiatively cool by hydrogen atomic-line excitation 
and efficiently form stars. The radiative transfer grid has 203^ cells. 

The H-ionizing photon luminosities per halo in our cosmic reionization simulations 
are assigned in the following way. Halo catalogues are discrete in time, because N-body 
density fields are stored every ~ 20 Myrs and the corresponding source (halo) catalogues 
are produced at the same time. A halo of mass M is assumed to have converted a mass 
M-{nb/D.m) ■ f* into stars, where /* is the star formation efficiencji||J. If each source forms 

II Not all the halos convert their mass into stars. Halos inside H II regions are assumed to be "failed" 
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Figure 5. Spatial slices of the ionized and neutral gas density from self-regulated 
simulation (IMSP; f2000.250S case), at redshifts z = 15.7 (top-left), 12.9 (top-right), 
9.9 (bottom- left), and 7.9 (bottom-right). The volume- weighted (mass- weighted) 
global ionized fractions at these redshifts are 6.4 x 10~^ (9.7 x 10~^), 0.12 (0.15), 
0.56 (0.62), and 0.99 (0.99), respectively. Shown are the density field (green) overlayed 
with the ionized fraction (red/orange), and cells containing active sources are shown as 
dots. Each slice has a thickness of 86 h^^ ckpc, while sources shown are from a thicker 
(~ 1.8/i~^cMpc) region. 



stars over a period of time At and each stellar baryon produces Ni ionizing photons per 

sources (i.e. their ability to form stars is suppressed by the photoionization which created the H II 
region) if their mass is below 10^ Mq. More realistically, source formation in halos of this mass 
range may have a gra dual dependence on halo mas s: as mass decreases, it is harder to form sources 
inside, and vice versa ( Mesinger and Diikstral . l2008h . Nevertherless, as halo population is dominated 
by lowest-mass halos, which are most vulnerable to photoheating in this mass range, we adopt our 
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stellar lifetime and is used only once per At, and if a fraction /esc of these photons 
escape into the IGM, the ionizing photon number luminosity of a halo of mass M will 
be given by 

Qi = T-. , (26) 

At ■ htuh 

where n is the mean molecular weight and rriH is the mass of a hydrogen atom. In 

this model, stars are produced in a burst, and they keep radiating with fixed Qi for 

At ~ 20 Myrs. It is noteworthy that the result does not depend on the detailed shape 

of the source spectrum, but only on a frequency- integrated parameter A^j. 

We calculate {L,y) for the LW background sources in a similar way. All H- ionizing 
sources also produce H2 dissociating photons, and their (L^) is also constant during the 
source lifetime At. In each succeeding time interval. At, new sources are identified with 
the halo catalogue for that time-slice in the N-body results and are assumed to emit 
radiation with constant (Lj^). This (Lj^) is proportional to {Qi/ f esc) {^un / ^i) ■, where 
the proportionality constant is the dimensional factor which indicates the frequency- 
integrated number of ergss"^ Hz~^ per LW photon released in the source spectrum. We 
then construct the future world lines of these sources across the source lifetimes. At. 
At a given observing redshift, we then draw past light cones from every grid point in 
the simulation box. When the past light cone of an observer intersects the world line 
of a certain source, we register the comoving distance Tos to that source and its flux 
contribution. 

We choose here a specific case from IMSP, the self-regulated reionization case 
"f2000_250S" with WMAP3 background cosmology. In this scenario, small-mass halos 
(10^ < M/Mq < 10^) host high-efSciency emitters with top-heavy initial mass function, 
or "IMF" (e.g. massive Pop III stars). The (hydrogen- ionizing) photon production 
efficiency, /^ = /*/csc^i, of these sources is approximated by /^ = 2000. On the other 
hand, large-mass halos {M/Mq > 10^) are assumed to host lower-efficiency emitters 
approximated by /^ = 250 (e.g. Pop II stars with Salpeter IMF). The simulation box 
has a volume (35 h^^ Mpc)^, with h = 0.73. The reader is referred to IMSP for more 
details. 

At this point, we need to deal with the fact that the simulation box (35 h~^ cMpc 
in this case) is smaller than the LW horizon (97.39Q;cMpc). We simply use a periodic 
box condition, attaching identical boxes around the domain of calculation. We use 5^ 
identical boxes in total, choosing the central one as the domain of computation. One 
may, instead, shift and rotate boxes before attaching them: however, this approach 
would not be able to remove the finite-box effect completely either. In the future, we 
will use a much larger simulation box, which will naturally reduce this effect. 

We parallelized our code using the message passing interface library (MPI) to 
calculate the evolution of the spatially-varying LW background on distributed-memory 
parallel computers. We used "Lonestar", a massively-parallel supercomputer at the 

simple prescription in this paper. 
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Texas Advanced Computing Center (TACC) at the University of Texas at Austin. 
It has a total of 5,200 cores of dual-core Intel Xeon 5100 processors and 11.6 TB of 
aggregate memory. The particular run of LW background calculation presented here 
- run separately from the cosmic reionization simulation - took about 15 hours of 
computing time on Lonestar, when we used 256 computing cores and about 1.5 GB 
memory per core. The numbers of halos on the radiative transfer grid of 203^ cells 
in the simulation box of 35 /i^^ cMpc were about 1.3 x 10'^, 7 x 10^, and 1.9 x 10^ at 
z ~ 15, 10, and 8, respectively. Note that the effective total number of sources for our 
LW background calculation is about 125 times the number of sources in the simulation 
box of this particular size, due to the length scale of tlw- 

3.2. Evolution of the Globally- Averaged Ionizing and Dissociating Radiation 
Backgrounds 

The growth and geometry of the ionized fraction of the universe during the EOR is 
illustrated by the selected time-slices shown in Figure [51 The corresponding evolution 
of the globally-averaged ionized fraction and the ionizing and dissociating radiation 
backgrounds is summarized in Figure O with several interesting features revealed. 
First, at 2 > 10, small-mass halos dominate the large-mass halos in contributing both 
ionizing and dissociating photons. This is easily understood in the framework of the 
standard ACDM cosmology, because the population of low-mass halos dominates over 
that of high-mass halos, both in numbers and in total mass. These low-mass sources, 
however, become "self-regulated" as the universe gets more ionized (since they are 
suppressed as sources in the H II regions) and are later almost fully suppressed at 
z < 10, while the collapsed fraction in the unsuppressible (higher-mass) sources only 
grows with time. Thus, both reionization and dissociation are dominated by high-mass 
halos at z < 10. Second, n-JriYSN becomes smaller at z < 10 than aX z > 10. This 
is simply due to the transition of major source type from Pop III to Pop II, because 

(A^i/iVLw)ni > (A^i/A^Lw)ii. 

We note that the reionization history depends on the adopted values of the 
efficiency parameters, /^, for low- mass and high- mass halos. The ionizing background 
is degenerate in /*, /esc and A^i, in fact, as long as their product /^ is fixed. This is not 
true for the dissociating background, however, which breaks the the degeneracy between 
/esc and Ai, in general. Because dissociating photons — with emitted energy ranging 
from 11.2 eV to 13.6 eV — are largely unattenuated by their own interstellar medium 
inside the source halos, their escape fraction from the source halo is essentially unity. 
Therefore, for a given /^, or a given reionization history, 

Jim, 21 oc 1/fesc- 

Hardness of the spectral energy distribution (SED) of the source, which can be 
characterized by A'i/A'LW; also affects Jlw, 21- For a given f-y, 

JUN, 21 OC (A'i/A'LW)"^ • 
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Assuming a top-heavy IMF, Pop III obje cts have (A^i/A^Lw)ni 



Tumlinson and ShuU 



2000 



15 (e.g. 



Bromm et al.l . 120011). Pop II objects with a Salpeter IMF, 
on the other hand, have (A^i/A^Lw)ii ~ 1 (jTumlinson and Shulll I2OOOI and references 
therein). 

The spatially-averaged mean intensity is aheady as high as (Jlw, 21) = 0.1 by 
z ~ 15, when the mean ionized fraction of the universe is only (x) ~ 0.02, if the 
fiducial value of /esc, in = 0.2 is used for Pop III objects. Jlw, 21 is not affected by /esc, 11 
until z ~ 15, because Pop II objects start to emerge only after z ^ 15. As described 
above, Jlw, 21 is proportional to 1/ fesc- For instance, if /esc,iii = 1 instead of 0.2, then 
(Jlw,2i) = 0.1 will be reached later after z ~ 15. (Jlw,2i) depends on both /esc,iii 
and /esc, II at 2; < 15, and if /esc,ii=l as well, then (Jlw, 21) = 0.1 at z ^ 13, or when 
(x) ^0.1. 
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Figure 6. Global history of cosmic reionization and the dissociation background: 
(Left panel) (top) Evolution of the photon to baryon ratios, where rijon is the total 
ionizing photon number density accumulated until given redshift z, while tilw is 
the "instantaneous" dissociating photon number density at z given by jilw = 
^/ii5cv hv'^'^ — l^ i2&cv ^'^^)- (bottom): Global mean ionized fraction {x) vs. 
redshift. (Right panel): Global mean dissociating intensity, {Ju^ 21)' multiplied by 
(/csc/0.2), vs. the global mean ionized fraction. 



It is useful to compare our numerical results for the space-averaged LW intensity 
with the homogeneous universe approximation. The global mean dissociating intensity 
(Jlw, 21) is easily predictable if the fraction of mass collapsed into stars and the number 
of dissociating photons created per stellar baryon are known. Consider a universe where 
sources are homogeneously distributed, with a collapsed fraction which is identical to 
that obtained numerically from our N-body simulation. The emission coefficient (in 
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erg s~^ Hz~^ sr~^ cm~'^) is then only a function of redshift as follows: 

JuAZs) = -r^^sPm{Zs)j^fco\\{Zs)f*, (27) 

An i Im 

where eu^{eTgs~^ Hz"^ g^^) is the emissivity at emitted frequency Ug, Pm{zs) is the mean 
mass density at redshift Zg, and /coul^^s) is the source halo collapsed fraction. Note that 
this is the proper emission coefficient: the comoving emission coefficient is given by 

= -;—^usPm,o-^^—fcoii{zs), (28) 

where Pm,o is the mean mass density at present. Finally, the (proper) mean intensity in 
this homogeneous universe is obtained by 

Jo ^ + Zs 

where the source redshift Zs is implicitly related to ros by equation ([8]). The emission 
coefficient j^, (-^s) is shown in Figure [TJ where the "stair-steps" reflect the fact that, 
as described in § 12.31 , we use halo catalogues that are discrete in time, which results 
in a discontinuous evolution of /cou as well. The resulting (Ji/)lw, homo, discrete plotted in 
Figure [8] agrees well with the globally-averaged value of {Ju)LW,sim from our simulations, 
as expected. For comparison, we also plot in Figure M the homogeneous approximation 
result when the discrete, time-stepped collapsed fractions, /cou, shown in Figure [7] are 
replaced by a smoothly varying mass function based on fitting the N-body results over 
time. The importance of the H Lyman line opacity in attenuating the LW photons is 
illustrated by the quantity ( J) lw, homo, thin, also plotted in Figure [HI the unattenuated 
mean intensity in the homogeneous approximation if we neglect H Lyman line opacity 
but take account of the ultimate horizon which corresponds to the distance from which 
Lyman limit photons at 13.6 eV redshift to the minimum energy (11.5 eV) of interest 
here for the LW dissociation rate. 

Which sources are the dominant contributors to J^, sources near to or sources 
far from the observer? Let us consider the same homogeneous universe as described 
above. If photons were not attenuated (/mod = 1), and the comoving emission coefficient 
jusi^s) remained constant over time, then spherical shells with identical proper thickness 
A ijT^) would contribute equally to J^,, as seen in equation (see [21]). This is not the 
case, however. The factor /mod increases as Tos decreases (Figure [3]), and /cou usually 
increases as Zg decreases as well. These two factors combine to make nearby sources more 
important. An opposite trend can occur, however, if we consider only the sources inside 
low-mass halos, because the evolution of their emission coefficient is not monotonic, as 
seen in Figure [71 That nonmonotonicity is not enough to completely offset the increase 
of /mod, however. Defining the fractional contribution from spheres of varying radius Tos 
as 

/(< ^os) = / L °^ j ju,{Zs) ■ fmod{ros) 
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Figure 7. Comoving emission coefficients j,^^ (zs) from Pop III objects (dotted) residing 
in low-mass (10^ < M/Mq < 10^) halos and Pop II objects (solid) residing in high-mass 
{M/Mq > 10^) halos. These are constant during the source lifetime At « 20Myrs, 
because halo collapsed fractions, obtained from the self-regulated simulation results 
(IMSP), are also assumed to be constant during At, which is also identical to the time 
interval between adjacent N-body outputs. 
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we have calculated /(< Tos) both for low-mass halo sources and high-mass halo sources. 
The results for the low-mass (suppressible) and high-mass (not suppressible) halos are 
plotted in Figure [91 In the case of Jlw.hm, for example, while all sources within tlw can 
contribute, ~ 80% of Jlw.hm is contributed by sources at r < [0.35 — 0.45]rLw In the 
case of Jlw,lm, the overall trend is similar to that for Jlw.hm, except at 2; ~ 8. Even at 
this redshift, however, an observer primarily needs to look back only to r ~ OAru^, since 
the high- mass halos dominate the total emissivity at late times, as seen in Figures [7] and 
[HI and dominate Jlw, as well, therefore. 

While our comparison of the simulated, globally-averaged LW background with 
that from the homogeneous approximation shows that the latter is good as long as 
we give it the correct, self- regulated, space-averaged mass function of source halos, 
only the simulations can tell us about the spatial variations in the LW background 
and their evolution, as well. We will show in the following section that a significant 
spatial fluctuation of the LW background does arise during the epoch of reionization. 
This, in principle, could induce a fluctuating feedback effect on the star formation in 
minihalos, thereby altering the apparent pattern of their clustering from that which 
arises gravitationally due to structure formation. 

Note that the possible effect by the finite size of the box, which is smaller than the 
LW horizon, does not seem to affect our quantitative conclusion too much. First, as 
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Figure 8. Evolution of the global average LW background intensity (in units 
of 10~^^ergcin~-^s~-^Hz~^ster~^) (J,y)LW,2i: {Ju)lw,21 (solid-cross), the intensity 
obtained by averaging the LW background field over the simulation volume at each 
redshift z, { J u)ysn, 21, ho-m (dot-dashed-square), the mean intensity in a homogeneous 
universe but with the emission coefficient depicted in Figure [7] (i.e. distrete in time). 
We also plot the average intensity, { J v) lw, 21, se-mi (long-dash), based upon our smooth 
fitting formulae (i.e. continuous in time) for the collapsed fractions of both low- 
mass (10^ < M/Mq < 10^) and high-mass {M/Mq > 10^) halos. For comparison, 
(>/i/)LH',2i,scmi, OT (short dash) is the intensity when we neglect the optical depth to 
H Lyman lines but take the horizon as the distance over which Lyman limit photons 
(13.6 eV) redshift to the minimum energy (11.5 eV) photons of our interest. 



described above, about 80% of the LW intensity conies from r < [0.35 — 0.45]rLW; which 
is just about the box size used. This imphes that the flucuation of LW intensity as well 
as its mean is contirubted mostly by nearby sources. Second, according to a suite of 
structure formation simulations we have performed, halo mass functions dn/dM from 
these simulations do not show too much suppression at a mass scale corresponding to 
the mass of the box, where most of the suppression of power is expected to occur. 
Mass functions from these simulations (in varying dynamical ranges), accordingly, 
connect smoothly when plotted on a single viewgraph. Because LW intensity originates 
from sources of reionization residing in cosmological halos, power of the LW intensity 
fluctuation at the box scale would not be suppressed just as the mass function dn/dM 
at the box scale is not. 



3.3. Spatial Fluctuations of the LW background 

Until now, the possible spatial fluctuations of the LW background have been neglected, 
due partly to the fact that tlw is a cosmologically large scale. One might naively 
have expected, therefore, that the LW intensity fluctuations inside the LW band photon 
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Figure 9. Fractional contribution to Ji, from sources inside a sphere of comoving 
radius ros (normalized to tlw) from the observer. This is separated into two categories 
by mass: {top-lejf) contribution to Jlw.hm, or the intensity Jlw only by sources in 
high-mass {M/Mq > 10^) halos. {bottom- left) contribution to Jlw.lm, or the intensity 
Jlw only by sources in low-mass (10^ < M/Mq < 10^) halos. {right) Fractional 
contribution to the toal Jlw by low-mass halo sources (LM, dot-dashed), high-mass 
halo sources (HM, dashed), and all sources (All, solid). Low-mass halo contribution, 
dominating Jlw at 2: ^ 20, becomes comparable to high-mass halo contribution at 
z ~ 14, and becomes negligible at 2; ~ 8. Considering both (All) contributions, about 
80% of Jlw comes from r < [0.35 — 0.45]rLw at all redshifts. This indicates that both 
isotropic and fluctuating components of the LW background are dominated by nearby 
sources at all times. 



horizon, or tlw; are negligibly small, since the fluctuations in the space-density of matter 
are small when averaged on such large scales at such an early epoch. Nevertheless, as 
we will show, we find huge fluctuations in the LW radiation field. How is this possible? 
In this section, we present our results for the fluctuations in the LW radiation field and 
explain their origin. 

In Figure [TOl we show the Jlw, 21 field on a planar slice (with fixed comoving 
coordinates) inside the simulation box at different redshifts. Jlw, 21 varies significantly 
in space at all redshifts. For instance, at z = 16.602, the overall variance of Jlw, 21 
is about two orders of magnitude. Figure [10] shows that as many as 3 contour levels 
are observed simultaneously in the same image plane at this redshift, corresponding 
to Jlw, 21 = 0.01, 0.1, and 1, respectively. This example clearly debunks any previous 
assumption that might have been made that only small fluctuations would exist in the 
LW background inside tlw- 

The PDF of Jlw, 21 values on the simulation grid, plotted in Figure [H], is also 
noteworthy. The volume- weighted distribution of Jlw, 21 is highly skewed. The deviation 
of Jlw, 21 from the global average, (Jlw, 21), is small when Jlw, 21 < (Jlw, 21) (or 6j < 0, 
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where 6j = (J — {J))/ {J))- The minimum fluctuation is (5j^mm ~ —0.5 at our starting 
redshift, approaching zero as time goes on. Roughly speaking, Jlw, 21 — {Jim, 21) when 
(5j < 0. In contrast, strong deviations are observed in regions with (5j > 0. Accordingly, 
most of the fluctuation involves Jlw, 21 > (Jlw, 21)- At all redshifts, Jlw, 21 shows a 
variation of about two orders of magnitude from the minimum to the maximum, and 
about an order of magnitude variation at the 99.73% level. 

The (sample) variance of Jlw, 21, o^^j = (^j), also plotted in Figure \T\\ does not 
evolve in a monotonic way. It starts out very large at 2; ~ 19.2 and slowly decreases in 
time until z k, 15, then increases again until z k. \2 (except for a brief, sudden drop 
at z ~ 13), and finally decreases again after z ^ 12. This limited range of redshift 
where a'^j reverses its evolutionary trend, 2; ^ [15 — 12], corresponds to the epoch when 
sources start to form inside massive halos with M > IO^Mq, which are not subject to 
suppression and "self-regulation" (see Figure [7j). We suspect that the small sudden drop 
of 0"} at z ~ 13, however, is just due to the inherent cosmic variance. 

The power spectrum of 6j, P{k), is shown in Figure [121 Large scales (small k) have 
more power than small scales (large k) do. The shape of P{k) is rather complex, and a 
simple power-law provides a poor fit. Nevertheless, if a local power-law fitting is used, its 
power index becomes steeper in time. The overall amplitude of P{k) decreases in time 
almost monotonically, except for an increase from z ^ 15 to z ^ 12, which is consistent 
with the trend seen in the evolution of aj in Figure fTTl As for the normalization constant 
of P{k), we follow the convention that the variance of 6j is given by 

^"^(^") = ^2l Pik)k'dk, (31) 

in the limit in which the size of the volume over which the average is calculated becomes 
infinitely large. 

What causes such huge fluctuations within the simulation box, even when its 
size is smaller than tlw? The answer is straightforward: radiation sources cluster on 
scales smaller than tlw, and their spatial clustering generates fluctuations in the LW 
background that are not washed out even after the fluxes from all sources within a 
distance tlw are added up. The original assumption that the LW background would 
evidence only small fluctuations misses this important ingredi ent. As simulated and 



noted by llliev. Mellema. Pen. Merz. Shapiro and Alvarezl (l2006l ) and others, patchiness 



of cosmic reionization, itself, strongly reflects the source-clustering effect. Sources 
cluster in high-density regions and will also produce a stronger LW background nearby, 
therefore. Such a correlation between the LW fluctuations and the matter density 
fluctuations is depicted in Figure [T3l 

Since reionization is "inside-out" according to these simulations (i.e. the high- 
density regions ionize first), there is also a correlation between the H II regions and the 
regions of higher-than- average LW intensity, as seen in Figure [TOl In fact, an animated 
sequence of maps like those in Figure [TO] shows that isocontours of Jlw,2i start out 
centered on the same density peaks where H II regions first appear. However, the 
isocontours of Jlw,2i expand more rapidly than the ionization fronts ("I-fronts") that 
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define the H II boundaries, overtaking the I-fronts and expanding beyond the H II 
regions. 

Finally, we focus our attention on the LW intensity field in the neutral regions 
alone. Some fraction of the Pop III stars, including the very first ones, are believed to 
have formed inside cosmological minihalos. Minihalos in ionized regions are less likely 
to have formed stars, however. If a minihalo formed in the ionized region, "Jeans-mass 
filtering" wou l d hay e meant that it formed of dark matter only, devoid of baryons (e.g. 



Shapiro et al.l . Il994l ). On the other hand, if a pre-existing minihalo found itself inside an 



ionized region that had expanded to overtake it before the minihalo ha d yet formed a 



star, such a minihalo w ould have photoevaporated in the ionized region (IShapiro et al. 



2004 



Iliev et al.l . I2OO5I ). Accordingly, star formation in minihalos in ionized regions was 
suppressed even more readily and further than in the low-mass atomic-cooling sources, 
halos with mass M > 10^. As far as star formation in minihalos is concerned, therefore, 
the neutral regions of the IGM are of particular interest. We find (see Figure [T^ that 
early-on, the PDF distribution of Jlw,2i in the neutral regions is very similar to that 
overall, except with the highest-flux tail of the distribution cut off, since those are the 
regions in the immediate vicinity of the ionizing sources, which are ionized first. The 
standard deviation of the PDF of Jlw, 21 decreases in time more rapidly for the neutral 
regions than for that overall. During the late stages of reionization, the values of Jlw, 21 
in these neutral cells converge to the average value, (Jlw, 21), and the fluctuations largely 
disappear. 

4. Conclusion and Future Prospects 

We have, for the first time, calculated the inhomogeneous background of II2 dissociating 
UV radiation caused by the same sources which reionized the Universe in a large-scale 
radiative transfer simulation of cosmic reionization. The UV continuum emitted below 
13.6 eV by each source was transferred through the IGM, attenuated by atomic H 
Lyman series resonance lines, to predict the evolution of the inhomogeneous radiation 
background in the energy range of the LW band of H2 between 11 and 13.6 eV. 
This required us to transfer the radiation from the many thousands of source galactic 
halos in our simulation, which resolved all halos of mass above IO^Mq in a comoving 
volume of (50Mpc)^, from each source to each of the millions of grid cells of our 
reionization simulation. To accomplish this, we developed a novel method to calculate 
the attenuation of LW band photons from individual sources by H Lyman series 
resonance lines, an otherwise prohibitively expensive, multi-frequency calculation, in 
a very fast way, instead, without an explicit multi-frequency operation. This was 
achieved by a grey opacity approximation, the "picket-fence" modulation factor, a simple 
function of the comoving distance between a source and an observer, which represents the 
frequency-averaged attenuation of LW band photons. We also explicitly accounted for 
the effect of the finite light-crossing time between sources and observers, by constructing 
the conformal space-time diagram of all reionization sources and observer grid cells and 
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finding the intersections between source world lines and observer past light cones. 

Our results here demonstrate that the rise of the cosmic LW background to levels 
which have previously been identified as the threshold for dissociating H2 and thereby 
suppressing star formation inside minihalos occurs well before the epoch of reionization 
is very advanced. Not only is this the case for the mean LW background intensity, as 
anticipated by earlier estimates based upon a homogeneous universe approximation, but 
it is even more the case for the inhomogeneous background. The first regions to form 
halos are the regions with the highest mean overdensity, and we show that this is where 
the LW background rises the fastest and is at the highest levels. This means that our 
assumption here, for simplicity, that reionization is dominated by the atomic cooling 
halos, and that minihalos are sterilized by the rising LW background before they can 
contribute significantly to the ionizing and, accordingly, the LW background, as well, is 
self-consistent. 

On the other hand, since there are also some minihalos that form far from the 
density peaks around which the halos cluster, we might also expect that there are some 
minihalos that form far from the peaks in the LW background, before their local LW 
intensity has risen to the threshold level for H2 suppression. In the future, it will be 
interesting to consider the possible role of these minihalos that form in places where the 
LW background is not high enough to suppress their star formation. 

Our result can be used for various applications. For example, fluctuating H2 
abundance in IGM and cosmological halos may be calculated inside our simulation box. 
This would at least require implementing reaction rates of neutral and ionic species of 
H, He, and H^ . One may have to run many small-box simulations similar to those of 



Yoshida et al.l ( l2003l ). in order to calculate H2 abundance and track source formation 
inside cosmological halos under fluctuating Jlw calculated in this paper. Similarly, once 
X-ray emitting sources are properly populated in the simulation box, a composite effect 
of negative (due to photodissociation by UV) and positive (due to partial ionization by 
X-ray) feedback may be studied as well. 

There is the possibility that source formation inside the more massive, atomic- 
cooling halos is affected by the LW background, too, because molecular cooling takes 
place inside these halos in the following way. When the g as cools from the ionized 



state in these halos, it does so out of equilibrium (e.g. IShapiro and Kang) 119871 : 



Kang and Shapiro! Il992l ): gas cools faster than it recombines, so even after atomic H 
cooling has reached its typical saturated phase at T < 10^ K, there still remains a 
significant trace amount of electrons. Gas-phase reactions can then create H2 with 
the help of these electrons, which can further cool the gas down to T ~ 200 K. 
Even though H-; at this stage can become self-shielded against the UV dissociating 



background ( Oh and HaimanI I2OO2I ). a strong enough background may nevertheless 
dissociate H2 and suppress star formation to some extent (HAR). The threshold LW 
intensity for such suppression inside these atomic cooling halos, (JLw)threshoid,atomic5 
is believed to be much larger than (Jlw) threshold for minihalos (HAR). This might, 
therefore, somewhat affect the history of cosmic reionization near the end of reionization. 
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when even (JLw)threshoid,atomic has been reached. 

Emissivity by the first stars that contribute to the cosmic near-infrared background 
(NIRB) may also be affected by the fluctuating LW background calculated here. 
Obser vations of a strong NIR B excess over the known foreground have been reported 



[e.g. iMatsumoto et al.l 120051). and this has been interpreted as a possible signature 
of the first stars (c.f. iMadaul |2006| and references therein), al though the true 
ident i ty of the contributing sources is currently under debate (e.g. iKashlinsky et al 



20071 ; [Thompson et al.l 120071 ). Our work will impact the interpretation of both the 



homogeneous and inhomogeneous components of the NIRB, because predictions of the 
first star formation should be strongly affected by our results. 

We argued in § [2] that we are justified in neglecting H2 self-shielding by the IGM 
in calculating the opacity to LW photons. This was justified by the fact that the 
concentration of H2 in the pre-reionization IGM was only ~ 10~^ and, as such, th2 ^ 1- 
However, we can now use hindsight to see that our neglect of th2 is even more self- 
consistent than that estimate would suggest. Hydrogen molecules in the IGM are 
quickly dissociated by the rising LW background. It is easily seen from Figures O 
[To] and [TT] that Jlw,2i ^ 10^^ at 2; ~ 20 in the vicinity of atomic cooling halos, and 
the mean value quickly rises to (Jlw,2i) — 10^^ and above. When the mean IGM 
prior to reionization (with temperature Tjqm — 369(1 + z^/lSb"^, hydrogen number 
density rin — 1-7 x 10"^ (1 + 2;)^/21^cm~^, and 2/H2 — 2 x 10"^ when there is no LW 
background) is exposed to the LW intensity Jlw,2i = 10^^ from z = 20 onward, for 
example, its molecular fraction drops to yn^ — 2 x 10^^ by ^ ~ 17. Ev en if H9. self- 
shield ing were marginally important with yn^ ^ 2 x 10^^ at high redshift (JRicotti et al. 



2OO2I ) . therefore, this will quickly become negligible when molecules are dissociated to the 
level ?/H2 — 2 X 10^^, which would occur at z ~ 17 from the LW background contributed 
by the atomic-cooling halo s alone. Peculiar mot io n of gas eleme i its wi ll weaken self- 
shielding even further (e.g. iMachacek et al.l 1200 ll ). iJohnson et al.l (120071 ) . on the other 
hand, claim that relic H II regions created by the first stars can recombine and generate 
abundant H2 before being exposed to other external sources, through the nonequilibrium 
H?-fo r mation mechari i sm de scribed above (IShapiro and Kang 119871 : iKang and Shapiro 



I992I ). Ijohnson et al.l (120071 ) find, however, that this boost of H2 abundance is delayed 
significantly when Ju/v, 21 ^ 10~^, which is easily satisfied in the mean IGM at z ~ 17 
and in the vicinity of atomic cooling halos at z ~ 20. This suggests that the proposed 
nonequilibrium enhancement of the H2 concentration in the IGM inside relic H II regions 
is not likely to provide a significant enough H2 opacity to shield the IGM through most 
of the EOR. 
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z=19.175 



z=16.602 




Figure 10. Simulation spatial images showing the isocontours of patchy reionization 
and the patchy H2 dissociating background on a planar slice through the box of volume 
(35/i~^cMpc)'^ at different epochs. The level of Jlw, 21 on the grid is depicted by 
different colors, with the range [10^'^ — 10^], shown on the inset of the top-left panel. 
On top of each Jlw, 21 color-map, contours of thick colored lines represent different 
Jlw, 21 levels (red, orange, blue, cyan, and green corresponding to Jlw, 21 =0.01, 0.1, 1, 
10, and 100, respectively). The black lines represent the ionization fronts characterized 
by a; = 0.5. 
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Figure 11. {left) Probability distribution function ("PDF") of Jlw,21 for the 203^ 
grid cells inside a box of comoving size 35 h~^ Mpc at different redshifts. Numbers on 
individual curves represent the corresponding redshifts. {right) The top panel shows 
deviation of Jlw. 21 from the globally-averaged (Jlw,2i)j expressed in terms of 1 + 5j. 
Around the curve of (Jlw, 21) (solid), contours containing 68.27% (dotted), 95.45% 
(short dashed), and 99.73% (long dashed) of the Jlw, 21 distribution are shown. The 
sample variance of Jlw, 21, cr,/ ^i is plotted together with the variance on the average 
, cr?j>, in the bottom panel. (T?j> is dominated by the Poisson error from the number 
of radiation sources rather than the number of simulation grids, because the former is 
found to be much smaller than the latter at all redshifts of our interest in our simulation 
box. 
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Figure 12. {left) Power spectrum P{k) of LW background fluctuations 5j at different 
redshifts. Numbers next to plotted curves represent corresponding redshifts. {right) 
P{k) at limited range of redshifts, z w 15 — 13.4. These show a reversed evolutionary 
track, increasing in time, compared to P(/c)'s plotted on the left panel, decreasing in 
time. Note that this trend can also be seen in Figure [TlJ Error bars on both plots 
represent variance of P{k) due to the finite number of wavenumbers and the finite 
number of radiation sources. At z = 19.175, there is only one radiation source in the 
box, and the corresponding power spectrum is roughly identical to the upper limit of 
P{k). 
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Figure 13. Correlation of Jlw, 21 and pm at different epochs, depicted by the number 
of cells at given 5j and (5p^, with uniform bin-size of log(l 



5p) and log(l + 5j^ 



Redshifts of these panels match those of Fig [TOl The inset in the top left panel shows 
the color scheme to depict the number of grid cells. 
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Figure 14. Probability distribution of Jlw. 21 in neutral cells only. Compare this plot 
to Fig [TlJ Decrease of overall area of PDF in time simply reflects the fact that the 
total volume of neutral regions decreases as cosmic reionization proceeds. 



